Josephson current in carbon nanotubes with spin-orbit interaction 
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We demonstrate that curvature-induced spin-orbit (SO) coupling induces a — n transition in the Josephson 
current through a carbon nanotube quantum dot coupled to superconducting leads. In the non-interacting regime, 
the transition can be tuned by applying parallel magnetic field near the critical field where orbital states become 
degenerate. Moreover, the interplay between charging and SO effects in the Coulomb Blockade and cotunneling 
regimes leads to a rich phase diagram with well-defined (analytical) boundaries in parameter space. Finally, the 
phase always prevails in the Kondo regime. Our calculations are relevant in view of recent experimental 
advances in transport through ultra-clean carbon nanotubes. 



The spectrum of quantum dots (QDs) defined in car- 
bon nanotubes (NTs) is four-fold degenerate owing to spin 
and valley symmetry. Recently, Kuemmeth et al l have 
demonstrated that the spin and valley degrees of freedom 
are coupled in NTs. This spin-orbit (SO) coupling breaks 
the four-fold degeneracy into two Kramers doublets (time- 
reversed electrons pairs). From a different perspective, 
NTs are interesting because they can support supercurrents 
when coupled to superconductors 2 5 . These supercurrents 
mainly result from resonant transmission through discrete 
states confined to the QD, the so-called Andreev bound 
states (ABS) corresponding to entangled time-reversed 
electron-hole Kramers pairs 6 . As both phenomena, SO and 
ABS, are related to time-reversed Kramers pairs, it is thus 
interesting to raise the following question: how are the 
ABS, and therefore the Josephson effect, affected by SO 
coupling in NTs? Here we address this question. Using 
various theoretical approaches we analyse this problem in 
all relevant transport regimes and demonstrate that the SO 
coupling is able to reverse the supercurrent, namely to in- 
duce a to 7r transition, even in the non-interacting regime. 

The valley isospin (r = ±) originates from the two 
equivalent dispersion cones (K and K') in graphene, aris- 
ing from time-inversion symmetry. When graphene is 
wrapped into a cylinder to create a NT, the valley degen- 
eracy leads to two degenerate clockwise and counterclock- 
wise electron orbits which encircle the NT. This degener- 
acy, together with spin, manifests in a four-fold shell struc- 
ture in the Coulomb Blockade regime 7,8 , as well as in a 
SU(4) Kondo effect in the strongly correlated regime 9 ' 10 . 
Furthermore, magnetic moments associated with these 
orbital persistent currents are remarkably large 11 which 
allows to perform detailed transport spectroscopy when 
an external magnetic field is applied parallel to the NT 
axis 91112 . The orbital motion of electrons also couples 
to a curvature-induced radial electric field. This creates 
an effective axial magnetic field Bso which polarizes the 
spins along the NT axis and favors parallel alignment of 
the spin and orbital magnetic momenta {K, f) and (K' , \) 
or antiparallel (K, l) and (K' , t) depending on the sign 




FIG. 1. (Color online), (a) Schematics of a NT coupled to super- 
conducting reservoirs. In the QD region, discrete Andreev levels 
form inside the BCS gap. The figure also show the K, K' orbits 
encircling the NT. (b) Energy spectrum of a NT QD for realis- 
tic experimental parameters (see supplementary info). All ener- 
gies are given in units of the BCS gap A = 0.25meV, such that 
Aso / A ~ 1.66) and referred to Ef which we take as the energy 
at which (K, f) and (K' , f) cross at B c » 0.52T (dashed verti- 
cal line), (c) Andreev bound states corresponding to the spectrum 
in (b). Black (orange) lines correspond to ABS calculated from 
the lowest (highest) Kramers doublet (each contributes with two, 
solid and dashed, ABS). (d) Critical current (units 2e A /h) versus 
gate voltage. The two peaks correspond to resonant Cooper pair 
tunneling through SO split Kramers pairs. 



of Ago- As a result, the fourfold degeneracy breaks into 
two Kramers doublets (time-reversed electrons pairs) sep- 
arated by an energy Ago 13 . Recent experiments 14 have 
shown this SO effect also appears in disordered NTs in the 
multielectron regime. 

The system we have in mind is shown in Fig. la. 
A QD NT with SO coupling is connected to supercon- 
ducting leads with BCS density of states. Owing to the 
superconducting pairing, electrons in the NT with en- 
ergies below the superconducting gap (A) are reflected 
as their time-reversed particle, a hole with opposite spin 
and momentum. This process, known as Andreev reflec- 
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FIG. 2. (Color online) Total (top), discrete (middle) and continuous (bottom) Josephson current (units 2eA/ft) as a function of phase 
and different B\\ (in Tesla) for Y = 0.1 A. At the highest magnetic filed the system has 7r -junction behavior, (b) Same as (a) near the 
0-7r transition at B\\ = B c = 0.52T. (c) ABS vs. 4> for different B\\. When B\\ > B c , the two inner ABS cross at E F — resulting 
in 7r behavior, (d) ABS versus V g for different B\\ =0, 0.5, 0.52 and 0.6 T, from left to right. At Bu = B c the two inner ABS are 
degenerate for all j V a \ < A. The n transition is robust as V g is varied (direction of the arrow) either above (e) or below (f) V a = 0. 



tion, leads to discrete states inside the gap, namely the 
ABS corresponding to entangled time-reversed electron- 
hole Kramers pairs. We model this system by an Anderson 
hamiltonian with s- wave superconducting reservoirs and 
with QD levels obtained from a NT model including SO. 
Green's function in Nambu representation are used to ob- 
tain the ABS and the two contributions to the Josephson 
current Ij = Ij 13 + Ij on of this model (full details are 



given in the supplementary info). The discrete part Iy s is 
due to Cooper pair tunneling through the ABS and can be 

written as if = f E Bl(2) /(^i(2)) with f(E) 
the Fermi-Dirac function. Namely, the derivative with re- 
spect to phase of the occupied ABS. The continuous part 
jcon j s fi ue to p ar ticle-hole excitations for energies larger 
than A. In the noninteracting case, the ABS can be ob- 
tained from 



(f f i TE m \ ( F +F | TE m \ r 2 AW(^/2) 

-^1(2) - £=Ft / ^1(2) + £ ±; H , X2 S^2 ~ U ' 
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where 4> is the phase difference between superconductors 
and r is the tunneling rate. The notation ia(2) indicates 
whether the Kramers doublet which contributes to the ABS 
is the ground (excited) state at B\\ =0 (Fig. lb). Impor- 
tantly, each Kramers doublet gives two solutions in Eq. (1) 
so in general we obtain four ABS. The two outer (inner) 
solutions correspond to Ei( 2 ) (Fig- 1°)- The results for the 
Josephson current are shown in Figs. 2a and 2b where a 
0-7r transition occurs for Bu > B c . The transition can be 
understood by studying the ABS spectrum as a function of 
4> for different B\\ (Fig. 2c). When B\\ > B c , the two inner 
ABS cross at Ep =0. Owing to this, the occupied ABS 
for Bn > B c belong to the same Kramers doublet (the 
one formed by (K, X) and (K' , I) which are, of course, no 
longer degenerate). Importantly, they carry supercurrents 
of opposite sign which leads to a negligible I^is- The main 
contribution is thus given by the continuum part which re- 
sults in 7r behavior 15 . In Fig. 2d we plot the ABS as a func- 
tion of gate voltage and different B\\ . At zero magnetic 



field, the SO splitted ABS show a diamond-like shape, 
similarly to spin-slit ABS due to Coulomb Blockade 16 ' 17 . 
As B\\ increases, the diamond closes and, ultimately, the 
two inner ABS become degenerate when B\\ — B c . When 
B\\ > B c , the ABS cross at Ep. After the crossing, the 
occupied ABS belong to the same Kramers doublet for a 
large range of \V g \ < A resulting in a n transition which 
is robust as the gate voltage is varied (Figs. 2e,2f). We 
include the effect of the Coulomb repulsion by first con- 
sidering the large gap limit, i.e., A — > oo, where the prob- 
lem can be mapped onto an effective low-energy model 
(U <C A) with a superconducting pair potential due to the 
proximity effect A^ = Tcos(<p/2). Direct diagonalization 
produces results for the ground state energy Ecs(<j>) an d 
trivially Ij = lj ls (in this limit, this is the only contribu- 
tion to the Josephson current). Owing to SO, the total spin 
S and orbital pseudospin T are no good quantum numbers. 
Instead, % d has a block diagonal form using the total pro- 
jections (S Z ,T Z ), with e T . s = e + 1/2tsAso, as a basis. 
For <j> = 7r, we find the analytical solution 
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- iA so < U < -e 
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A so + 3U for -\e - \A so <U<-\e- \A so , (S Z ,T Z ) 
fmU<-U-±Aso,(S z ,T x ) = (0,0) 



(0,0),(1,0),(0,1) 
= (±1/2, ±1/2) 



(2) 



The ground state for arbitrary has to be calculated nu- 
merically (Fig. 3a shows the phase diagram for 4> = 0). 
Nevertheless, it can be shown (by comparing with the ap- 
proximate boundaries obtained by perturbation theory in 
A D , red lines in Fig. 3a) that for large U the ground state is 
always (S Z ,T Z ) = (±1/2,^1/2) with energy E GS {<p) = 
2e - l/2( A /4r 2 coa 2 (0/2) + (A so + 2e + 3U) 2 + 3U). 
While we cannot identify this state with a n phase, it is 
likely that the inclusion of quantum fluctuations, by con- 
sidering a finite gap, will stabilize the system towards 
this phase. Indeed, cotunneling corrections for (T <C 
A), present tt phases. This can be shown by employ- 
ing second-order perturbation theory in T (namely fourth- 
order cotunneling processes, see supplementary info) 1 . In 
this limit, we find a supercurrent Ij = I c sin(^) such that 
the overall sign of I c governs the or 7r-character. In par- 
ticular, the 0-7r -junction transition takes place at the value 
of e corresponding to the resonant condition e — A so /2 = 
Ep = 0, with a tt phase for e < A so /2, such that the tran- 
sition can be tuned by a gate voltage. Numerical results are 
shown in Fig. 3b. 

Beyond cotunneling, higher order tunneling events lead 
to Kondo physics. Here, we consider the large-?/ limit 
(supplementary info) where simultaneous fluctuations in 
the spin and orbital quantum numbers lead to a highly 
symmetric SU(4) Kondo effect (for a Kondo temperature 



T K ,su(i) » A )- When T K . SU {i) > A so , we find 

jdis = eA sin(0) 

J " 2H ^=± [(1 + ^) 2 + !][(! + 7?a) 2 + COS 2 (|)]' 

(3) 

with a = 2T ^° (4) ■ When T KiSU{4) < A so , only the 
lower dot level participates in producing an SU(2) Kondo 
state. In the limit T KiS u(2) > A, the ABS are simply 
Ei = ±Acos(0/2), namely the ABS of a single contact 
with unitary transmission. The corresponding supercurrent 
is If s = sinO/2), with |0| < tt 19 . Fig. 3c summa- 
rizes these results. For both symmetries, the Josephson 
current always exhibits a 0-junction behavior but the mag- 
nitude strongly depends on Aso, as shown in Fig. 3d. For 
Aso = 0, we recover the results of Ref. 20 . 

In closing, we have demonstrated that SO coupling in- 
duces a — 7r transition in the Josephson current through 
a QD NT coupled to superconducting leads. Our calcu- 
lations, which cover all relevant transport regimes, non- 
interacting, Coulomb Blockade, cotunneling and Kondo, 
determine in a precise manner the conditions for the tran- 
sition in terms of system parameters which can be tuned 
experimentally. Our predictions are relevant in view of 
recent experimental advances in transport through ultra- 
clean NTs with SO coupling 1 . Furthermore, most of the 
physics discussed here is inherent to the rich behavior that 
ABS show in the presence of SO coupling. We therefore 
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FIG. 3. (Color online), a) Phase diagram in the large gap limit and = 0. For large enough U, the ground state is (S Z ,T Z ) = 
(±1/2, Tl/2) (dark pink region), b) Critical current versus level position in the cotunneling limit. The critical current undergoes 



a 0-7T transition when e 



Aso/2 = 0. c) Discrete Josephson current (in units of ^) versus 
Josephson current for 4> = it versus SO coupling. As the system changes from SU(4) to SU(2) Kondo symmetries, Ij l 
to maximum. 



in the Kondo limit, d) Discrete 
goes from zero 



expect that tunneling spectroscopy of individual ABS, like 
in the experiments of Ref. 17 , may also reveal the effects de- 
scribed here. Microwave spectroscopy of excited ABS 21 is 
one further experimental example where our findings may 
be tested. 
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Appendix A: NT Model 

We consider a single wall NT whose low energies can 
be described by expanding the momentum near the Dirac 
points of graphene H — hvp(k y T 3 ® o\ + k x r ® <t 2 ), 
here vp is the Fermi velocity, r 3 is a Pauli matrix acting on 
isospin (K, K') space (with eigenvalues r = ±1) whereas 
the Pauli matrices <j\ and 02 act in sublattice space (the 
two carbon atoms in the primitive unit cell of the graphene 
honeycomb lattice). k x and k y are the momenta along the 
NT axis and circumferential direction, respectively. The 



eigenvalues of Ho are E (k x ,k y ) — ±hv F ^Jk 2 + k 2 . Im- 
posing periodic boundary conditions, k y is quantized as 
k y — 2tv/3D (lowest mode), where D is the NT diame- 
ter and v depends on the type of tube. In the following, we 
will consider small bandgap tubes parametrized as k y — 
rk g . We also include a magnetic field Bn applied paral- 
lel to the NT axis. B\\ induces an Aharonov-Bohm flux 
$ AB = B\\irD 2 such that k y = rk g + $ab/D§ , with 
^0 = h/e being the flux quantum. Besides this orbital 
shift, B\ 1 also induces the standard Zeeman shift in the spin 

sector Hz = hgfiB^\\ T o ® °o ® s ^ with s 3 being a Pauli 
matrix (eigenvalues s = ±1) describing the spin projection 
along the tube axis. Finally, the SO coupling term has the 
form Uso = {^so T 3 ® °i ® s 3 + A^ g t 3 ® a ® s 3 ), 
which includes off-diagonal A^ and diagonal A° so 22,23 
terms in sublattice space . The eigenvalues of the full H = 

Hq+Hz+Hso read E s>T (k x , k y ) = ±hv F ^Jk 2 x + k^ + 
s(TAg C + \g[j,BB\\), here Ag has been absorbed in k y 
as k y = rk g + &ab/D&o + sAg /hvp. Finite interval- 
ley scattering Ak ,k 1 introduces anticrossings in the spec- 
trum when spin polarized orbital states are degenerate (not 
shown). 
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1. QD Bound states 



The total (low-energy) Hamiltonian for a quantum dot 
carbon nanotube with spin-orbit coupling can be written 



H(x) = hv F (k y T3®(Ji + k x TQ®(j2)®s + (&iT3 (8) (7i (g) s 3 + A r 3 ® a ® s 3 ) + ^.9/iB£||To®cro<g>s 3 + V"(af) (Al) 
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FIG. 4. Energy diagram of a CNT-QD with electrostatic gates 
inducing a potential V(x) along the nanotube axis. Bound state 
energy E (green dashed line) obeys a criterion E g < Et < Vq + 
E a . 



Here Vf is the Fermi velocity, t 3 is a Pauli matrix acting on 
isospin (K, K') space (with eigenvalues r = ±1) whereas 
the Pauli matrices <j\ and a 2 act in sublattice space (they 
account for the two carbon atoms in the primitive unit cell 
of the honeycomb lattice describing graphene). k x and k y 
are the momenta along the NT axis and circumferential 
direction, respectively. The term V(x) describes the po- 
tential induced by the electrostatic gates and is defined as 



a simple step potential of the form 

V(x) 



Vq, \x\>1 
0, \x\<l 



We use the following ansatz 24 ' 25 for the electronic wave- 
functions ip T (t) defined in different intervals (see Fig. 1) 



ViiAx) = Ce- q " x rKM*,c/v 



with momenta k x = 



f^-) - k 2 
\hv F J V> 



Tk„ 



^± (for the lowest, m = 0, subband). 



The corre- 



sponding energies are e T . s = ± l ^jf L ^(Lk x ) 2 + {Lk y ) 2 + 



(A r + \g^ B B\{) s, E 
and zl , , = ± rk y~ ik * 

k x ,k y ,c/v y/k^+kl 



. Hvp 



A t + \gn B B\\) s 
with L = 11. Here, the 



subscripts c and v correspond to conduction and valence 
bands, respectively. The energy levels e T S are found from 
the continuity of the wavefunction ip T (x) at all potential 
steps. That is, 



Ae -q*e^k v -iq a ,c/v^ =Bie -ik x e^k v ,k*,c/v^j +B2e ik x e^l y ,-k x ,c/v^j ^_ if,J(-£) = ip}j{-£) (A3a) 
Ce~ q * 1 ^k y -^,c/v^j = Bl& ik x l ^k y ,k x ,c/v^ + B2e -ik x i ^k y ,-k^c/v^j ^_ il) T n {£) = ipj n {£) (A3b) 



Using Eqs. (A3a), we solve for B\ and B2 and get 

(A4) 



B x \ ±iAe ( e %,- iqx , e /v e z k y ~k x ,c/v 



*/ LK X j * / fC x -T K y \ k y ,-iq x ,c/v k y ,k x ,c/v 

where +(— ) belongs to the conduction (valence) band. Also, from Eq. (A3b) we have 



Z k V ,k X ,c/V + ^ Z ky,-k X ,c/v 



^k y .iq x , c /v Bie ik " i + Bie~ ih * 1 
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Substituting Eq. (A4) into Eq. (A5) yields 



Tk y +q x y/(kl-q*Hk*+kl 



-- (rky sin 2k x £ — k x cos 2k x t) — sin 2k x £ 



\f k y $ 



(rky sin 2k x l + k x cos 2k x t) 



■■ sin2fc T i' 



(A6) 



Simplifying the above equation gives 



sin 2k 7 



y/{kl-ql){kl + kl) 



(k y sin 2k x t + k x q x cos 2k x t) 



(A7) 



which yields 



tan2fc T ^ 



k x q x 



k x q x 



^(k 2 y -q 2 x )(k 2 x + k 2 y )-k 2 y \E-V \\E\/{hv F f-kl 



(A8) 



Fig. lb in the main text shows a typical energy spectrum for realistic NT parameters. In particular, we use SO coupling 
values A = 0.26meV and Ai = 0.053meV such that Aso — 0.4156meV. Rest of parameters: k g = — 0.09nm _1 , 
L = 100 nm, V = 3.95meV, hv F = 526.57meV nm . 



2. Including superconducting leads 

The QD CNT coupled to superconducting leads is mod- 
elled as an Anderson Hamiltonian coupled to s-wave su- 
perconductors with BCS density of states. This Hamil- 
tonian can be written in second quantization as H = 
He + Wd + Wt where 



a—L I R,k,r,s 



-Eh* 



e^cL^c^., +h.c. 



(A9a) 



U D = y~] £ TS d\ s d TS + U ^2 n TS n T > sl (A9b) 

T,S (T,S)^(T' ,S') 

H T = ]T ( V ^ikr S drs + h. C ) , 



(A9c) 



a—L/R : k,T : s 



where c' akrs creates an electron on lead a e {L, R} with 
energy and with quantum numbers k, s, and t corre- 
sponding to the wave-vector, spin and orbital degree of 
freedom, respectively. A is the superconducting gap and 
(f> = — 4>ji is the phase difference. d TS is the op- 
erator that annihilates an electron on the dot with energy 
e TS (where the dependence on gate voltage is implicit). 
U denotes the intra- and inter-orbital charging energy and 
n TS = d\. s d TS represents the occupation operator for the 



dot levels. The last term describes tunneling by means 
of energy-independent tunneling amplitudes V a leading to 
tunneling rates T a — ir\V a \ 2 p (p is the contact density of 
states). 

Appendix B: Calculation of Andreev bound states and the 
Josephson current by using the Green's functions technique 

A powerful technique to obtain the total Josephson cur- 
rent through the system described above is the Green's 
function method where all physical quantities can be writ- 
ten in terms of the Green's functions 

g r d a (t,t>) = ((d,S)) r ' a = TiQ(±tTt') ([d(t),£(t')} + 



g<(t,t') = ((d,S))< =i(S(t')d(t) 

Q>{t,t') = {{l$))> = -i(d(t)d\t') 



(Bla) 



Owing to the superconducting pairing, these Green's func- 
tions are matrices containing anomalous components. In 
the following, we write these matrices using the following 
the Nambu bispinors 



Cak+t 

C 'k 
ak— 4. 

C f - 



and d ■■ 



(B2) 
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Using the standard relation z(( A, B) ) + ( ( [H , A] , B)) = (\A,B] + ), together with the following commutation relations 



[H,d TS ] = -e TS d TS -Udrs^ nT ' s ' ~ VgCgki 

t' ,s' a.k 

[H,dU] = Sf S d\ s + Ud) fS Y^ n f's' + VaC lkfs 



a.k 



[H, c akTS ] = -ikCakrs + sgn(s)Ajf a cj^__ - V a d TS 
[H, c\- k ] = ^ k c\, - S gn(s)A-^c akTS + V a d\s, 



we obtain 



/d +t n +t \ 



n-id 



l u -i 



d-^ri-t 



g r a ^((c ak ,d^y = a 3 v a ((d,d^y 

where n TS = J2' T > s > n T's' (the prime in the summation means (t', s') ^ (r, s)) and 



z + £-i 



Z — £_f 



z + £+u 



and 



v a = 



vj 



( z - £ fe A Q 
V 



(B3a) 
(B3b) 

(B3c) 
(B3d) 



(B4a) 
(B4b) 



(B5) 



/I 



0-3 



(B6) 



-1/ 



Owing to the presence of Coulomb Interactions, U ^ 0, 
the equations of motion for the Green's functions, Eqs. 
(B4), cannot be closed and we need some approximations 
which we discuss next. 



Green's functions and self-energies can be obtained. In 
particular, the retarded ones read: 



1. Non-interacting limit (7 = 

a. Retarded Green' s functions 



((d,S)y= fgj-'-s; 



(B7) 



In the noninteracting case, U = 0, the equations of mo- 
tions can be closed such that analytical expressions for the 



where 



^ r = ^2va 3 g r ak a 3 v 

a,k 



f3 d 

(3 Q cos(0/2) 



[3 cos(0/2) 



[3 d /3 o cos(0/2) 
/3 cos(<£/2) [3 d J 



(B8) 



withT = 2-k Pn {Q)V 2 and 



_M_ 



^/uj 2 -AI 
A Q 

■ sgn(^) A a 
^ 2 -A2 



if M < A 
if M > A 
if |w| < A 
if U > A 



Using At, — Ar = A, V L = V R = V, and ct> L 
<f>/2, the retarded Green's function reads: 



(B9a) 



(B9b) 



where 



r>4 



z + e^ + Tpa -r/3 cos(<^/2)\ 
-r^cos^/2) z-e +1 - + T/3 d ) 







j_ / * + e +i + Tf3 d -Tp cos(0/2)\ 
D- \ -T/3 cos(0/2) z - e_ t + J / 



= e+t + r A*)( z + + r /3d) - ( r & cos(0/2)) 2 
D- = (* - e-t + rfl,)(* + e +i + Tfa) - (T(3 cos(0/2)) 2 . 



(BIO) 



(Blla) 
(Bllb) 



b. Andreev Bound States 



When \u>\ < A, the Andreev bound states can be determined from the poles of the Green's function. Namely, we just 
need to solve the determinant equation DetlQ^u;) -1 ] = 0. Using Eq. (B7) we obtain the following equation: 



= D+D- = 



Explicitly, 



(w - e +t + Tf3 d )(uj + £ _ 4 + Tf3 d ) - (Tf3 cos(0/2)) 2 (w - e_ t + rfl,)(w + e +i + r/3 d ) - (r& cos(0/2)) 2 



(B12) 

0. 

(B13) 

At this point it is important to note the full equivalence of the Green's function method with the Bogoliubov-DeGennes 
Hamiltonian method (indeed, the Green's function has precisely Bogoliubov-DeGennes structure). The Andreev bound 
states give rise to delta-function contributions in the spectral density. The weights can be found by the residues of the 
Green's function at these poles. Explicitly, 



11 E 2 

((d,S))]^ = -nJ2Qb+S(uJ - E 2 ) 



.33 



.34 



Si 



-TT^Q b _S(0J - 
Ex 



where 



Z\y± = lim 

ui— >-E 2 (i) 



0,22(44) 



( W ) 



J 0, 12(34) 



lim 



14 M 
( W ) 



and 



#2(1) - £±t + 



(i) 



TA 2 



#2(1) + St-1 + 



FE 2 (i 



(i) 



V /A2 ^( 1 ), 



(A 2 -w 2 )VA 2 



2w - e ±T + e Ti 



2Tcj 



VA 2 - uj 2 



r 2 A 2 cos 2 (0/2) _ 
A 2 - K 2 _ 

^ ^2(1) 

2wr 2 A 2 cos 2 (0/2) 
(A 2 - w 2 ) 2 ' 



(B14a) 
(B14b) 
(B14c) 
(B14d) 

(B15a) 
(B15b) 

(B16a) 
(B16b) 



Eq. (B16a) corresponds to Eq. (1) in the main text. 



c. Josephson Current 



ing to the Josephson effect, this expression contains a 



jdc 



The current through a given lead a can be written as 



e d ^t^ with n o 



E 



Ow- 



9 



dissipationless component (nonzero current at zero bias 
voltage) when there is a superconducting phase differ- 
ence. Thus, the Josephson current can be extracted from 
the general by just studying the limit of zero ap- 
plied bias voltage I a = 



\v dc =o 



2e» 

h Jl 



Tr 



d{n a ) 
~ C dt 



V da =0 



|v dc =o- Using the 



tion methods, one finds that the Josephson current can be 
expressed as 



(B17) 



nonequilibrium Green's function and the equation of mo- where 



( Pd 

fa-*** 

V 



Pd 



Pd Poe 



(B18a) 



P e 



= E v^i^v = re(M - A)/( W ) 

fe 

withT = 2 7 rp A r(0)V r2 and 



V 



if |w| < A 
if M > A 



,+»</ 



/ 



(B18b) 



ifM<A 

ifM>A. 



(B19a) 



(B19b) 



One important advantage of this method is that the Josephson current can be easily split into two parts I = his + Icon- 
The first part is the so-called discrete contribution and corresponds to the Josephson current carried by Andreev Bound 
states. The second term, the so-called continuous part I con , corresponds to the current given by the continuous spectrum 
of states above the gap. Both expressions can be written analytically as: 



hi, 



2 ' r J ge(A-H)/H 7 =^= sin(^,/2)a \((d,&)) r 21 + «d,c?»; 2 + ({d,$)y 43 



Icon 



h J 2tt 

2eT f duj 

IT 



^e(H-A)/( w 



VA 2 - oj 2 

sgn(w)A 
Vw 2 - A 2 



sin 



43 



where, again, A^ = Ajj = A, Vj, = Vr = V, and 4>l = —(j)R = <j)/2. Explicitly, 

.f(£ 2 )A 2 , x . /(-Ei) A 2 



e r 2 

his = --z-sin(i^) 



E 



(A 2 - E 2 )D' + (E 2 ) 



E 



7rfi 



sin(0) / 9(M - A) 



/HA 2 , 
(w 2 - A 2 )' 



(A 2 - Ef)D'_(Ei) 
1 1 



£>+(w) D-(w) 



«d>^»34j 

(B20a) 
(B20b) 



(B21a) 
(B21b) 



After some algebra, the discrete contribution can be rewritten as 

2e ^ <9£i( 2 ) 



(B22) 



-El(2) 



which is the expression discussed in the main text. 
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2. Cotunneling regime 



Expressions in the cotunneling regime can be obtained by lowest (second order) perturbation theory in T 18 . Starting 
from the expression for the current 



fe,r,s 



where 



(B23) 



(B24) 



we perform a standard thermodynamic perturbation expansion in the tunneling and obtain the Josephson current in the 
lowest non-vanishing order (fourth order in Ht) as 



h 



1 



rP 



•0 



dn dr 2 dT 3 (T T (H T (T 1 )nT(T2)n T (T 3 )nT a )) l 



3! Jo Jo Jo 

The Josephson current must involve two and two H^, which can be chosen in three ways, and hence 

h 



(B25) 



f dn j dr 2 ( dr 3 (T T (H+ a (r 1 )H+ gi (T 2 )H T jT 3 )H Ta )) 
o Jo Jo 



(B26) 



where we have used that in order to have Cooper pair tunneling, the H J must belong to the opposite junction. Next, if we 
choose the valley and spin of the last as, say, (+, t), it then means that the other / K T carries (—,4-)- m the same way, 
the valley and spin of the two can be chosen. All in all, we thus obtain 



f dn f dr 2 [ P dr 3 (T T {n+ afl Ari)n+ &T ,Ar2)n^ s (r3)n TaTS )) 

Jo Jo Jo 



(B27) 



At arbitrary B\\, the Josephson current can be written as I a = I c sm(<p), where the critical current reads 

eT 2 



where 



and 



N(e) 



,00 ,00 A A 

= -/ dE dE' 

J a J a y/E 2 - A 2 y/E' 2 - A 2 



C(E,E',e,s) 



(B28) 



(B29) 



C(E,E',e,e) = - 



rP rP 

. v / dn / dr 2 \ dT 3 T L {E,r 3 )T R {E\n-T 2 )Q{T 2 -T 3 )Q{r 3 -nY r ^ e{r ^ Tl) 

' e TS Jo JO Jo 

(B30) 



with e 



TS C TS- 



The functions T a are related to the the anomalous Green's and given by 
functions of the leads, which are defined as 

^ ak+ (r,r') = - (t t cI- (r)cl k+t (r')) (B31a) ^(r,/) - ^-(r,r') = ^ ak (r,r') 

A a e 



^ ak -(T,T>) = (Trcl k+i (r)ci k _ t (r')\ (B31b) 



2E, 



ak 



-F a {E ak ,T-T r ) (B32a) 
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where Here, first performing the imaginary time integration and 

then taking the approximation exp[— /3E] w 0, the func- 

T a {E ak ,T) = e - E " kM -2 cosh (E ak T)n F (E ak ) tion C(E, E', e, e) is given by 

(B33) 

Throughout, we assume low temperatures such that 
Al/rP 3> 1, and we thus approximate 



(B34) 



C(E,E',s,e) = 



-f3e 



+ 



(E + E')(E + e)(E> + e) (E + E')(E + e){E' + e) 



/3(e-e) 



+ 



(E + E' — e + e)(E — e)(E' - e) (E + E' + e - e){E - e)(E' - e) 



. (B35) 



In general, the integrals in Eq. (B29) have to be evaluated well below the Fermi level. Then, at T — only the lowest 
numerically. For example, let us assume that all levels are level contributes so that the critical current is given by 



eT 2 
2hir 2 



[N(s +l ) + N(e- t )} 



.^M 2 (-e +i /A) forB||>0 
■S E M 2 (-6_ t /A) forB|,<0 



where 



M 2 (e) 



I ""L 



dv 



vV - 1 vV - 1 (u + v + e + e){u + e)(v + e) 



(B36a) 



(B37) 



Even in this simple case, the integral cannot be solved analytically. This is in contrast with the limit B\\ = 0, where 
further analytical progress can be made. Assuming, for simplicity, = S-_i = e d + Ago/2 and = e + i = 
e d — Ago/2, the or tt character of the junction can be extracted by the overall sign of the critical current which reads 

Ic=& [N(s +t ) + AT( £ _ t )] with 



N(s +t ) = { 



2M(e +t /A)/A fore d >+A so /2 

for -A so /2 < e d < +A so /2 

for e d < -A S o/2 



and 



2M(e_ t /A)/A for e d > +A so /2 
-M(-e_ t /A)/A for -A so /2 < e d < +A so /2 
-M(-e_ t /A)/A for e d < -A so /2 



Here, the dimensionless function M(x) defined as (x > — 1) 

/>oo /*00 ^ 

/ du I dv — 



M(x) 



Vu 2 - 1 Vv 2 ^! (u + v)(u + x)(v + x) 



(B38a) 



(B38b) 



(B39) 



can be expressed by 



M(x) = 



(tt/2) 2 (1-x) 



arccos" x 



where the analytic continuations of arccos x = i ln(x + 
Vx 2 — 1) for x > 1 is understood. The function M{x) 



x(\ — x 2 ) 



with x > — 1 
(B40) 
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is always positive, it diverges at x — » — 1 + , and then 
smoothly decays for increasing x with the asymptote 
A4(x) <~ n 2 /4:X 2 for x — > oo, which allows to extract an- 
alytical boundaries for the — tt transition. In particular, 
we can establish the following criteria 

1. If both levels and £_+ are above the Fermi level, 
the Josephson current is positive, i.e., 0-junction. 

2. If both levels are below the Fermi level, it shows a 
tt -junction behavior. 

3. If the level e+f is above the Fermi level and the other 
level e_| is below the Fermi level, it is again a ir- 
j unction. 

For Ago = 0, the function N(ed) can be written as 

'2M(e_ T /A)/A fove d >0 



N(e d ) 



-i.M(-e_ t /A)/A fove d <0 



(B41) 



an the critical current is given by I c — j^*r&N(ed) = 
^ [46(e d ) - ®{-e d )\ M(\e d \/£S), in agreement with 
Ref. 20 . 



3. Kondo regime 



We study the Kondo regime in the large-L/ limit by 
means of the slave boson method. Using this language, 
Eq. (A9) can be rewritten as 



n SB = Yl Zk4 tkTa c ak T. - E (^^IkrAm + hM ) 

a—L/ R.k.r.s 



a,k,r 



+ J2e TS fUrs + ^= E (yJ akTS tfUs + h.c) +Kr£flf TS + tfb-l \ , (B42) 

t,s V a=L/R,k,T,s \t,s ) 



where the physical fermionic operator is written as d\ s = 
fl s b. The pseudofermion operator /t g creates a state with 
spin s and isospin r and the slave boson operator b annihi- 
lates an empty state. It can be shown that this mapping is 
exact provided that the constraint 



J2fLfrs+b^b=l 



(B43) 



In order to obtain self-consistent equations, we replace 
(6+) by y/Nb* and obtain 



a.k 



A\bf 



(B45) 



is fulfilled (A is a Lagrange multiplier which enforces this 
constraint). Note that the hybridization element is rescaled 
into V a -> V a /\fN. From the equation of motion for the 
boson field b\ we have 



a,k,r,s 

(B44) 



Eq. (B45) constitutes a set of self -consistent equations to- 
gether with the constraint 



jj(fHt)*3f(t)) + \b\ 2 = jj (B46) 
This mean field approximation, which neglects charge 
fluctuations, is reliable in the deep Kondo regime where 
only spin fluctuations are relevant. 



In the frequency space, the equations become 

du , 



a.k 



2i:i 
1 

N 



-Tr 



V a a 3 ((c ak ,p))< 
du> 



A- 



2m 



-Tr 



0"3 



r 
r 



i 

N 



(B47a) 
(B47b) 
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where 



(V a b 



V a 



V a b 



V n b 



(B48) 



V 



V a bJ 



and T = r 



. At this point, we have to calculate the lesser Green's functions. To do that, we note that the mean-field 



Hamiltonian is given by 



cx—L I R,k,r,s 



a.k.r 



+ J2^s.fLfrs+ ( V «b*ci kT J TS + h.c)+k(\b\ 2 -l) (B49) 

T,s a=L/R,k,T,s 



where e TS = e TS + A. The retarded Green's functions can be then written as 



«/,/ t » r = 



gf 



where 



f z — e 



+t 



— r, — 1 

8/ = 



z + e-i 



Z — £_f 



(B50) 



(B51) 



and 



z r = J2 va ^3V = -r 



p d /3 o cos(0/2) 
/3 cos(</>/2) (3 d 



\ 



[3 d /3 o cos(0/2) 
/3 o cos(0/2) fa J 



The lesser Green's function is given by 

«/,/ t » < = «/,/ t » r S < «/,/ t » a = -/M (((/,/ t )) r - ((/,/ t )) a ) 
Eqs. (B47) can be further simplified as 



-Tr 



tw< 



A ./ 2ttz 

1 f du! 



N / 2tr 



-Tr 



«/./ t »; 



0"3 



r 

A-=0 
f _ 1 

+ r ~~ a' 



(B52) 

(B53) 

(B54a) 
(B54b) 



Using these equations we obtain analytical expressions for the renormalized parameters e and T, from which we can 
extract the Kondo temperature and the position of the Kondo resonance. The slave boson mean field hamiltonian in Eq. 
(B49) is quadratic such that we can use the techniques in the previous sections to obtain the Andreev bound states and the 
Jospehson current. In what follows, we discuss these quantities in different regimes. 

First, let us consider the deep Kondo regime in the absence of the spin-orbit coupling A S o = 0. Then, the effective 
level is given by e = T K SU ^y Using this fact, the Andreev bound states can be written as 



E b /A = ±4 



li 2 +T 2 cos 2 ((/>/2) 



e 2 + Y 2 
Employing Eq. (B22) we find 



± 



\ 



T£,su(4) + T lsu { ,) cos 2 (g/2) h 

1 K,SU{A) + 1 K,SU(A) V Z 



eA 



sin(0) 



eA >— ^ sin(0) 
* s ~ 2^2fi ^1 + cos 2 (0/2) ~ V^fr v 71 + co S 2(^72) 



(B55) 



(B56) 
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On the contrary, for SU(2) model the effective level is e = Then, 

E b /A = ±cos(0/2) 

which implies 

eA 

his = ~y sin (0/ 2 ) 



(B57) 



(B58) 



Second, we consider what happens in the presence of the spin-orbit coupling. In this case, for T K> su{i) ^ Ago, the 
effective levels become s + f = S-i = T K< su(A) + Aso/2, £-f = £+4. = T K> su(i) — Ago/2, such that 



E b+ /A = ± 



\ 



P +t + T 2 cos 2 (0/2) 



? +t + f 2 



( 1+ 2T^) 2 +^ 2 (^) 

\ \ Z1 K,SU( 



E b -/A = ± 



\ 



? t + r 2 cos 2 (</>/2) 

?2 



+ r 2 



= ± 



\ 



f 1 A so ^ _i_ l 

\ 2T Jf,SC7(4) / 



(B59a) 



(B59b) 



which yield 



_ eA ^ 



sin(0) 



2ft 



(B60) 



r,=± 



( 1 + ^2T^) 2 + 1 ( 1 + '/57^) a + 00^(^/2) 



Next, we study the opposite limit T K ^su(i) "C Ago- ln this case, only the lower level participates in the Kondo physics. 
Thus, 



£-+ = £+i = for T K , SU(2 ) > A 
In this case, the Andreev bound states become E b -/A = ± cos(0/2) such that 

eA 

his = -^-sin(>/2), 



(B61) 



(B62) 



With 101 < 7T 19 . 
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